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Gutzwiller study of extended Hubbard models with fixed boson densities 
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We studied all possible ground states, including supersolid (SS) phases and phase separations of 
hard-core- and soft-core-extended Bose-Hubbard models with fixed boson densities by using the 
Gutzwiller variational wave function and the linear programming method. We found that the phase 
diagram of the soft-core model depends strongly on its transfer integral. Furthermore, for a large 
transfer integral, we showed that an SS phase can be the ground state even below or at half filling 
against the phase separation. We also found that the density difference between nearest-neighbor 
sites, which indicates the density order of the SS phase, depends strongly on the boson density and 
transfer integral. 

PACS numbers: 03.75.Hh, 05.30.Jp, 05.30.Rt 



I. INTRODUCTION 

The supersolid (SS) phase that simultaneously exhibits 
a superfluid (SF) phase and a density-wave order has 
been studied, since it was first proposed theoretically [l|- 
3]. Although its existence is still controversial, a non- 
classical moment of inertia of solid 4 He in Vycor glass 
that suggests its existence was recently reported |4{. 

SS phases may also exist in optical lattices. Cold 
atoms in optical lattices are new experimental systems in 
which theoretically modeled Hamiltonians can be more 
distinctly simulated experimentally with no disorders. 
For instance, an optical lattice system, as described by a 
Bose-Hubbard model 0, 0] , distinctly showed the phase 
transition between the SF phase and the Mott insulator 
(MI) phase Q. Moreover, recent observations of the long- 
range dipole-dipole interaction in 52 Cr atoms [§[ may 
lead to the realization of the SS phase. 

The simplest lattice model with a long-range interac- 
tion is the extended Bose-Hubbard model with a nearest- 
neighbor interaction: 
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H 



H kin + H int + H int , 
(i,3) 



u 



y]ni(rH - 1), 



(1) 



Here t, U, and V denote the transfer integral between 
nearest-neighbor sites, the repulsive intra-site interac- 
tion, and the repulsive inter-site interaction between 
nearest-neighbor sites, respectively. Furthermore, a, (aj) 
is the annihilation (creation) operator at the site i. In 
this study, we assume that the lattice is bipartite and 
consists of sublattices A and B and that the number of 
nearest-neighbor sites is z. The extended Bose-Hubbard 
model that prohibits (allows) multiple boson occupations 
at one site is called the hard-core (soft-core)-extended 



Hubbard model. 

Several analytical and numerical methods exist for 
studying the ground state of the extended Bose-Hubbard 
model. The strong-coupling perturbation theory [9( has 
been applied to obtain the phase boundary between the 
SF and non-SF phases. At least in the absence of the 
nearest- neighbor interaction [13, [ll|, the phase bound- 
ary determined by this theory agrees perfectly with that 
determined by quantum Monte Carlo (QMC) simulations 
[l2l - [T7l ]. However, this theory can neither distinguish be- 
tween the SF and SS phases nor describe discontinuous 
transitions even if they exist. In contrast, the mean field 
(MF) ap pro ximation and equivalent Gutzwiller approxi- 
mations [ll, [l^] can distinguish between the SF and SS 
phases and describe discontinuous transitions [20l - |22j |. al- 
though their validity is limited for high dimensions. 

For the hard-core model at half filling (N = 1/2), 
both the MF approximation j23rl25j for the mapped XXZ 
model and two-dimensional (2D) QMC studies [Ifl HH 
identified the first-order SF-solid transition and no SS 
phase. The MF approximation [23[ also showed that the 
transition from the SF phase to the checker-board solid 
phase is of the first order, and at the transition point 
fi = //*, the SS phase; solid phase; and PS phase in 
which the SS phase and the solid phase coexistx are all 
degenerate. 

On the other hand, later QMC studies [13, [H| showed 
that the SS phase is unstable against PS away from half 
filling. 

For the soft-core model, both Gutzwiller approxima- 
tions [Isl, and QMC simulations [HI, H(| have shown 
that checker-board SS phases can be stable against the 
PS. However, the details of the two types of study differ. 
The Gutzwiller studies showed stable SS phases for both 
broad filling and half filling, but the PS issue is not a 
concern in these cases. In contrast, the 2D QMC studies 
did not show the SS phase at half filling but identified SS 
phases above half filling; however, the system energy as a 
function of boson density was found to be concave below 
half filling, which indicates a possible PS instability pal ]. 
A later 2D QMC study [3l| showed more distinctly that 
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the SS phase is unstable against PS below half filling by 
exhibiting the negative compressibility k — dN/dfi < 0. 
This situation is somewhat similar in one dimension (ID). 
Both QMC [32I. [33j and density-matrix renormalization 
group [34H36I ] studies have been performed in ID. The 
SS phase was found only above half filling as in the 2D 
case but PS did not occur [H, HH . 

In a three-dimensional (3D) cubic lattice, Yamamoto 
and coworkers (37| found some region of the parameter 
sets where the SS phase is stable below half filling in 
QMC calculations, which is consistent with a MF phase 
diagram on the N-t/U plane for V = U/6. Hence, it 
is interesting to study the 3D system further. On the 
other hand, other previous calculations assumed a grand 
canonical ensemble and did not directly include the possi- 
bility that the system can be separated into two phases. 
As a result, the phases comprising the separate phases 
were not directly shown. However, if we start with a 
canonical ensemble and explicitly include the possibility 
that the system can be separated into phases, we can 
precisely describe the phase diagram with a fixed boson 
number for the entire system. 

Such calculations might also be easily considered us- 
ing the grand canonical ensemble if we explicitly as- 
sume that the total free energy of the system is Etat = 
7 -Ess + (1 — l)E so \. Here 7 is the ratio of the areas of 
the SS phases, and we neglect the surface energy between 
the SS and solid phases. However, if we minimize the free 
energy at the chemical potential /z, we can obtain only 
7 = or 1, but not an intermediate value < 7 < 1 
except for fj, = fi* at which the SS and solid phases 
are degenerate [Essi^*) = -^soi (/■**)]• That is, the solid 
and SS phases cannot be in the ground state at the same 
time, and a separate phase including two states cannot 
be obtained for /i 7^ /x*. To obtain a separate phase, 
we must tune the chemical potential from /1 to /1*; fur- 
thermore, we must try to find a 7 that satisfies the boson 
number condition 7A*ss + (1 — "f)N so i = N. These cal- 
culations become more complicated when the number of 
possible phases is large (SF, SS, MI, and several solids). 
This issue, in which PS can be obtained only on the phase 
boundary point(s) or curve(s) in the phase diagram but 
not for a finite region, generally occurs when we employ 
intensive variable(s) such as chemical potential in ther- 
modynamics. The most popular case is the gas-liquid 
transition, where PS occurs only on the phase boundary 
curve in the pressure-volume plane for a given temper- 
ature. However, if we employ particle number, volume, 
and internal energy as thermodynamic variables, which 
are extensive quantities, we can obtain a finite region of 
PSs on the particle density-energy density plane. 

In this study, we study the ground-state properties of 
the extended Bose-Hubbard model on a bipartite lattice 
based on the Gutzwiller approximation in a canonical 
ensemble; hence, we do not have to tune the chemical 
potential of PS. We study a wide parameter regime and 
present three main figures for the phase diagram on the 
particle density-nearest-neighbor interaction plane. Our 



phase diagrams for small, intermediate, and large trans- 
fer integrals differ significantly. The other main topic 
is the solid order parameter 8n — \Na — Nb\ [N^my. 
the boson density on the A(B) sublattice], which is also 
of interest and also depends strongly on the transfer in- 
tegral. We employed the linear programming method 
[38j to minimize the total energy for a particular boson 
number. Following this method, we can simultaneously 
determine the PS region, the phases that comprise the 
separate phase, and the ratio of each phase to the entire 
system. 

This paper is organized as follows, Section II intro- 
duces our calculation method, and Sees. Ill and IV 
present the hard-core and soft-core Bose-Hubbard mod- 
els, respectively. Section V describes the effect of an 
improved calculation on the energy of the solid and MI 
phases. Section VI discusses the conclusions based on 
our results. Finally, the appendices explain the details of 
several perturbative calculations, the results of which are 
compared with the numerical results from the Gutzwiller 
approximation. 



II. CALCULATION METHOD 



We employ the following Gutzwiller approximation for 
the SS and SF phases: \& = Y\ t $i. Here, $i is a varia- 
tional wave function at the site i. We assume a bipartite 
lattice with sublattices A and B and a checker-board 
symmetry for the SS phase. The variational wave func- 
tion is assumed to be $i = ^a(b) if the site * be- 
longs to the A(B) sublattice. $a(b) 1S written as a lin- 
ear combination of states \n) with n bosons as — 
J2 n c A(B)n\ n )- The variational parameters ca(b)u are de- 
termined so as to minimize the energy expectation value. 
We use Powell's method [38] to numerically optimize the 
variational parameters. If the result of the optimization 
shows CAn 7^ CBn and ca(b)u ^ &i{j)n i the phase is SS; 
however, if CAn = CBn 7^ §im the phase is SF (i and 
j are non- negative and positive integers, respectively). 
The first inequality for the SS phase corresponds to the 
existence of a finite density difference between sublattices 
A and B Sn — \Na — Nb\ (i.e., the checker-board density 
order). Here Na(b) = n \ c A(B)n I is the expectation 
value of the boson density at the A(B) sublattice. The 
second inequality implies that the phase is not solid. For 
the soft-core model, in principle, the sum of n must be 
taken from n = to 00; however, we take the sum from 
n = to n = 9, which is sufficient for our calculation 
when N < 1. For the hard-core model, by definition, we 
consider only the two states \n = 0) and \n = I). 

The energy expectation value per site is calculated by 
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the Gutzwiller variational wave function as 



III. HARD-CORE MODEL 



E 



(H) 

(H kin ) + (Hg t ) + « t ) 
—zt E v 7 "- + 1 



A.B n 



V 



E E^- 1 )!' 



A.B n 



I 2 + ^N A N B .{2) 



Furthermore, we assume that the energy per site of the 
MI phase with N bosons per site is 

Em = ^N(N-l) + ^N 2 = ^ foriV = 1(3) 

and that of the solid phase with Na {Nb) bosons per site 
on the A(B) sublattice is 

^ s = T E ^(N—^ + ^-NaNb, (4) 

t = A.B 1 



where N a{b) is a non- negative integer. In particular, for 
the solid phase when Na — 1 and Nb = (the Si 
phase), Eg 1 = 0; for the solid phase when Na = 2 
and N B = (the S 2 phase), E S2 = U/2. These en- 
ergies correspond to those obtained from the Gutzwiller 
variational wave function with CAn — CBn — 5nN 
for the MI phase and CAn = SnN A and cbu — S n N B 
(Na 7^ Nb) for the solid phase. Because we neglect the 
surface energy between different phases by assuming the 
thermodynamic limit, the system's total energy per site 
is 



E u 



for the boson number condition 

Not = Yj^ iNi - 



(5) 



(6) 



Here {-EJ and {JVJ (where % = SF, SS, MI, and solids) 
are the energies and boson number densities, respectively, 
of all possible phases and {7^} represent the of area or 
volume ratio of the phase i in the entire system. Follow- 
ing the linear programming method [HI, we can mini- 
mize E to t as a function of {7^} |39(. Only one 7^ value 
(corresponding to the uniform phase) or two 7^ values 
(corresponding to PS) are automatically chosen to be 
nonzero as a result of minimizing E to t because only one 
additional condition, the boson number condition, exists. 
For simplicity, however, we neglect the possibility of the 
separated phase that splits into the SF and SS phases, 
which is highly unlikely. Hereafter, we represent the sep- 
arated phase consisting of phases X and Y as the PS(X 
+ Y) phase. 



In the hard-core-extended Hubbard model, we prohibit 
multiple occupations at a site and define the Gutzwiller 
variational wave function as 



A(B 



) = C A(B )o|0) + C4( B )i|l). 



(7) 



At half filling, |c j4 ( jB ) | 2 + |c^_(_b)i| 2 = 1 according to 
the normalization condition of the wave function, and 
\cai\ 2 + \cbi\ 2 — 1 according to the boson number 
condition. If we set x — cao, then cai — cbo = 
\/l — x 1 and cbi = x. Hence, the system energy per 
site is 



E = 



zV 

-ztCAQCAlCBQCBl + — \cai\ 2 \cbi\ 2 



zt 



zV 



x 2 



(8) 



If V > 2t, then x 2 = or 1 and the phase is Si; 
however, if V < 2t, then x 2 — \ and the phase is SF. 
Hence, the SF-Si phase transition occurs at V = 2t 
and is discontinuous. This result agrees with those of 
previous studies (iTJ, |4l[ . 

Away from half filling, we can obtain the SF-SS phase 
boundary by a perturbative calculation setting ca(b)u = 
c n + 5cA(B)m where c„ is the optimized value when we 
assume the SF phase (cau — CBn) and 5ca(b)tl repre- 
sents infinitesimal quantities that describe the possible 
SS instability (see Appendix A.l for details). The result 



is 



t 



N 2 + (1 - N) 2 
N(l - N) 



(9) 



That is, the energy of the SS phase is lower (higher) 
than that of the SF phase for V > V c [V <V C ). Vg s is 
invariant under N O 1 — N because of the particle-hole 
symmetry of the hard-core model. We also obtain the 
phase boundary between the SF phase and the PS(SF + 
Si) phase, which consists of the SF and the Si phases by 
another perturbative calculation, assuming that the ratio 
of the Si phase to the entire system is infinitesimal (see 
Appendix A. 2 for details). The result agrees with Eq. [9] 

To determine the ground state for V > Vc, we nu- 
merically compared the energy of the SS and PS(SF + 
Si) phases by using the Gutzwiller variational wave func- 
tion as explained in the Introduction. We found that the 
PS(SF + Si) and SS phases are degenerate for V > Vc 
within the numerical error. The phase diagram is shown 
in Fig. 1. The numerical phase boundary agrees with Eq. 
[9] within the numerical error and the phase transitions 
are continuous as the perturbative calculation assumed. 
These results are consistent with those in Ref. 23] which 
showed that the SF-Si transition is of the first order and 
the SS, Si, and PS(SF + Si) phases are all degenerate 
at the transition point fi = fi*. 

We also confirmed that both 7s 1 (the ratio of the Si 
phase) in the PS(SF + Si) phase and Sn = \N A -N B \ in 
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FIG. 1: Phase diagram of the hard-core model. In the PS(SF 
+ Si)/SS phase, the SS and PS(SF + Si) phases are degener- 
ate. The two-dot-dashed line at half filling indicates the solid 
phase; the dashed line at unit filling indicates the MI phase. 

the SS phase arc finite for V > V c [Figs. 2(a) and 2(b)]. 
As N approaches 1/2, jg 1 increases rapidly as a function 
of (V— Vc)/Vc, because at half filling, no SS phase occurs 
and the SF-Si phase transition is discontinuous. 



IV. SOFT-CORE MODEL 

In this section, we study the soft-core-extended Hub- 
bard model. We employ the Gutzwiller variational wave 
function and optimize its variational parameters numer- 
ically. Hereafter, we call this the full numerical calcula- 
tion^); when we do not refer to the calculation method, 
the result was obtained by the full numerical calculations. 
We also perform perturbative calculations that limit the 
Hilbert space of the Gutzwiller variational wave function 
and include partial numerical calculations. We compare 
these perturbative calculations with the full numerical 
calculations. In subsection A, we examine the SS phase 
between the SF and solid phases at half filling. In sub- 
section B, we examine the case away from half filling, 
which is the main part of this paper. Here we show that 
both the phase diagram on the zV/U — N plane and 
Sn = \Na — Nb\ change qualitatively from small zt/U 
to large zt/U. 

A. HALF FILLING 

In this subsection, we primarily examine the SS phase 
at half filling. We compare the critical value of V for 
the SS-Si transition with that of the SF-SS transition 
because if the former is larger than the latter, we can 
obtain the SS phase between the SF and Si phases. 

By a perturbative calculation with an infinitesimal SS 
component added to the Si phase, we obtain the criti- 




(V-V c )/V c 

(b) 

FIG. 2: (a) V dependence of ys 1 [the ratio of the Si phase 
in the PS (SF + Si) phase] for zt = 1. (b) Sn = \N A - N B \ 
in the SS phase for zt = 1. Vc is the critical value for 
the SF-PS(SF + Si) or SF-SS transition at each value of N 
(zVc = 4.250, 2.762, 2.167, 2.002 for \N — 1/2| = 0.3, 0.2, 
0.1, 0.01, respectively). 

cal value of V for the SS-Si transition by a power-law 
expansion of zt/U: 

zV^ = 2zt + 2 Z -± + 2 Z -±+0(^L) (10) 

(see Appendix B.l for details). We can also obtain the 
critical value of V for the SF-SS transition 

by another perturbative calculation, in which an infinites- 
imal Sn = \Na — Nb\ is added to the SF phase (see 
Appendix B.2 for details). Because Vg F ~ ss < Vg s ~ Sl , 
the SS phase is possible for intermediate V satisfying 

ySF-SS < y < ySS-S^ 
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FIG. 3: Phase diagram of the soft-core model at half fill- 
ing. Solid curves represent numerical SF-SS and SS-Si phase 
boundaries. Dashed (dot-dashed) curve represents the SF-SS 
(SS-Si) phase boundary obtained by a perturbative calcula- 
tion using Eq. HIlfEq. fit))) . 




FIG. 4: V dependence of the density difference between 
sublattices A and B Sn = \Na — Nb\ in the SS phase of 
the soft-core model at half filling. V c is the critical value for 
the SF-SS transition at each value of zt/U (zVc/U = 0.219, 
0.469, 0.730 for zt/U = 0.1, 0.2, 0.3, respectively). 



We verified these results numerically. Figure 3 shows 
the phase boundaries obtained by the above-mentioned 
perturbative calculations (dashed and dot-dashed curves) 
and the full numerical calculations (solid curves). A sta- 
ble SS phase appears between the SF and Si phases. 
Both the SS-Si and SF-SS phase boundaries determined 
by perturbative calculations agree well with those deter- 
mined by the full numerical calculations (especially at 
small zt/U as expected). Figure 4 shows the interaction 
dependence of Sn = \N A -N B \. Here, zVc = zV^ F ~ ss 
is the critical value for the SF-SS transition at each value 
of zt/U. A finite Sn less than unity shows the density or- 
der of the SS phase, whereas Sn — 0(1) corresponds 
to the SF(Si) phase. We found that Sn continuously be- 
comes finite at the SF-SS phase transition. It changes 
more rapidly for smaller zt/U. This demonstrates that 
the SS phase disappears and the SF-Si discontinuous 
transition occurs in the hard core limit zt/U — > 0. 

On the other hand, 2D QMC calculations at half filling 
[25l I2H ] showed that the SF phase directly transitions to 
the Si phase and no SS phase appears. However, zt/U 
might not be sufficiently large to allow the SS phase to be 
easily found there, and a QMC simulation (as a matter 
of course, not only 2D but also 3D) for a large zt/U has 
a possibility for finding the SS phase. 

B. AWAY FROM HALF FILLING 

In this section, for comparison with the results of the 
full numerical calculations, we calculate the results of 
the following three perturbative equations for the phase 
boundary (see Appendices A, B.2, and B.3 for details). 
Perturbation 1 is Eq. [51 which assumes a limited Hilbert 
space with \N = 0) and N = 1) and is the same as 



that used in Sec. Ill; perturbation 2 is Eq. IB161 which 
yields the SF-SS phase boundary; and perturbation 3 is 
Eq. |B22 which yields the SF-PS(SF + Si) or SF-PS(SF 
+ S2) phase boundary. For perturbations 2 and 3, we 
employ a limited Hilbert space with | N = i), where i — 
0, 1, and 2. Hence, results obtained by perturbations 2 
and 3 are expected to be better than those obtained by 
perturbation 1. However, perturbations 2 and 3 require 
a numerical calculation to minimize the SF energy in the 
limited Hilbert space. 

We begin by analyzing the phase diagram for a small 
transfer integral zt/U. Figure 5(a) shows the phase dia- 
gram for zt/U = 0.03, which resembles that of the hard- 
core model. Namely, the upper bound of zV/U for the 
SF phase has a minimum value near half filling, which is 
similar to that of the hard-core model with particle-hole 
symmetry. PS occurs for large zV/U below half filling, 
and the separated phase is the PS(SF + Si) phase, as 
in the hard-core model. However, the SS phase appears 
above half filling. These properties agree with those ob- 
tained in a 2D QMC study with zt/U = 0.08 [Hj]. 

Figure 5(a) plots the result of perturbation 1 (dot- 
dashed curve) . It agrees almost perfectly with the full nu- 
merical calculation (solid curve) except at large N where 
the SF-SS phase boundary disappears because the PS(SS 
+ S2) phase appears there [see also Fig. 5(b), which ex- 
pands Fig. 5(a) around zVj U = 1]. Both perturbation 2 
for the SF-SS boundary curve and perturbation 3 for the 
SF-PS(SF + Si) boundary curve yield almost the same 
results as perturbation 1. Therefore, they also agree al- 
most exactly with the full numerical calculations except 
at large N. 

Figure 5(b) shows that PSs also appear above half fill- 
ing but the region is very small compared to that ob- 
tained by 2D QMC The appearance of the PS(SS 
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FIG. 5: Phase diagram of the soft-core model for zt/U = 
0.03. (a) Entire phase diagram. Solid curve represents the 
SF-PS(SF + Si) or SF-SS phase boundary. Dot-dashed curve 
for the SF-SS phase boundary obtained by the perturbative 
calculation in Eq. [9] cannot be distinguished from the solid 
curve except at large N, where the SF-SS transition curve 
disappears, (b) Expansion of (a) around zV/U = 1 and 
N = 1. Solid (dashed) curves for the phase boundaries 
show continuous (discontinuous) phase transitions. In both 
(a) and (b), the two-dot-dashed lines at half filling and unit 
filling indicate the solid phases (Si and S2, respectively), and 
the long-dashed line at unit filling indicates the MI phase. 



+ S2) phase and the intricate structure of the phase di- 
agram are non-trivial; however, the two SS phases that 
sandwich the PS(SS + S2) phase have different charac- 
teristics: one having a smaller zV/U resembles the SF 
phase and the other having a larger zV/U resembles the 
S2 phase, as we will see below through Sn = \Na — Nb\. 
Namely, SF, SS, S2, and the PS comprising these phases 
are almost degenerate for zV/U ~ 1 above half filling. 

Figure 6 confirms that the SS phase overcomes the 
PS(SF + Si) phase above half filling because the critical 




0.02 0.04 

zt/U 

FIG. 6: Critical nearest-neighbor interaction Vc for the SF- 
SS transition (solid curves) and the virtual SF-PS(SF + Si) 
transition (dashed curves) above half filling. Note that both 
solid and dashed curves appear for N = 0.6, although they 
are very close. 



value of the nearest-neighbor interaction Vc for the SF- 
SS transition is always smaller than that for the virtual 
SF-PS(SF + Si) transition. Here, the virtual SF-PS(SF 
+ Si) transition was obtained by setting the variational 
parameter of the Gutzwiller variational wave function to 
CAn — c-Bn- The difference between the curves for the 
SF-SS transition and the virtual SF-PS(SF + Si) tran- 
sition disappears in the hard-core limit of zt/U — > 0, 
where the SS and PS(SF + Si) phases are degenerate, as 
discussed in Sec. III. 

Figure 7(a) shows the V dependence of Sn obtained by 
the full numerical calculation. Figure 7(b) shows an ex- 
pansion of the region around 1.001 < zV/U < 1.002. 
Note that Sn has a discontinuity at zVq/U ~ 1.001, 
because the discontinuous SS-PS(SS + S2) transition 
occurs there. Interestingly, the SS phase exhibits a 
small Sn for zV/U < 1, whereas it exhibits a large 
Sn ~ 2N (N A ~ 2N and N B ~ 0) for zV/U > 1.002 
after a rapid increase in Sn in a narrow PS(SS + S2) 
region (1.001 < zV/U < 1.002), where the curves of 
Sn for different N are almost indistinguishable [sec Fig. 
7(b)]. Because Sn ~ 2N, Sn is larger for larger N at 
zV/U > 1.002. In contrast, Sn is smaller for larger N for 
zV/U < 1, demonstrating that the critical value of zV/U 
for the SF-SS transition is an increasing function of N 
above half filling. The fact that Sn ~ 2N for zV/U > 1 
indicates that the SS phase is similar to the S2 phase, in 
which N A = 2 and N B = 0. The fact that Sn ~ for 
zV/U < 1 indicates that the SS phase resembles the SF 
phase. Note that zV/U — 1 is also the phase transition 
point for t — at unit filling: the MI phase (i.e., no- 
density-order phase) for zV/U < 1 becomes the S2 phase 
(density order phase) for zV/U > 1. That is, even for a 
finite but small t and N < 1, the characteristics of the 
ground state change rapidly near zV/U ~ 1. 
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FIG. 7: V dependence of Sn — \Na — Nb\ in the SS phase of 
the soft-core model for zt/U — 0.03. Sn is larger at smaller 
N for zt/U < 1.000, discontinuous at zVc/U ~ 1.001, and 
smaller at smaller N for zt/U > 1.002. (b) Expansion around 
the PS(SS + S) phase (1.001 < zV c /U < 1.002). InthePS(SS 
+ S) phase, Sn is a sharply increasing function of zV/U, and 
the curves of Sn for different N cannot be distinguished. Thus, 
the curved arrows describe the increase in Sn at each value of 
N when we enlarge zV/U. 



Next, we study the case of an intermediate zt/U value 
of 0.3. Figures 8(a) and 8(b) show the phase diagram; 
the PS region above half filling exists, but again is very 
small. The phase diagram clearly departs from that of 
the hard-core model with the particle-hole symmetry. In- 
terestingly, for an intermediate zV/U, the SS phase over- 
comes the PS even below half filling. This agrees with a 
recent 3D QMC simulation with zt/U = 0.33 [37]. Fig- 
ure 9 shows Sn ( 7Sl ) below half filling in the SS (PS(SF 
+ Si)) phase. The SS-PS(SF + Si) transition is discon- 
tinuous: for N — 0.3, 0.35, and 0.4, Sn (7sJ is finite 
(zero) in the SS phase and discontinuously becomes zero 
(finite) in the PS(SF + Si) phase. For N = 0.25, Sn is 



zero in the entire figure because the SS phase does not 
exist. 

Returning to Fig. 8(a), we see that the results from 
perturbative calculations (dot-dashed curves) agree well 
with those from the full numerical calculations (solid 
curves). Assuming the Si phase for the solid, perturba- 
tion 3 is in almost perfect agreement with the SF-PS(SF 
+ Si) phase boundary below half filling. Perturbation 
2 is also in excellent agreement with the SF-SS phase 
boundary except for the region of large N ~ 1 where 
the SF-SS transition becomes discontinuous, and thus 
cannot be described by the perturbative calculation. 

Figure 8(b) shows an enlargement of Fig. 8(a) around 
zV/U = 1.1 and N ~ 1, where a small PS(SF + S 2 ) 
phase appears. The curve of perturbation 3 in Fig. 8(b), 
which is calculated assuming the S2 phase for the solid, 
seems to be far from the results obtained by the full nu- 
merical calculations; however, the difference is in fact 
approximately 1% at most. 

Figure 10 shows the V dependence of Sn above half 
filling for zt/U = 0.3. As in Fig. 7 for zt/U = 0.03, 
Sn ~ 2N is larger for larger N at zV/U > 1.1, whereas 
Sn is smaller for larger N for zV/U < 1.1. Hence, as 
in the case of small zt/U, the characteristics of the SS 
change drastically from SF-like to solid (S2)-like around 
zV/U ~ 1 when we increase zV/U . Note that in contrast 
to the case of zt/U = 0.03, Sn shows no discontinuity 
except at N = 0.9 because the SF-SS transition is con- 
tinuous. As a result, the two Sn curves for two different 
N values intersect smoothly around zV/U ~ 1.1 except 
at N = 0.9. For N = 0.9, Sn is zero for zt/U < 1.105 
(because the phase is SF, at which Sn — 0) and discon- 
tinuously becomes finite for zt/U > 1.105, because the 
SF-SS transition is discontinuous there. 

Finally, we examine a large zt/U( = 1). Figure 11 
shows the phase diagram. Only the SF and SS phases 
exist, and there are no PSs. The SF-SS transition is con- 
tinuous throughout the figure. The critical zV/U value 
between these phases is a decreasing function of the bo- 
son density N. This may be explained as follows. For 
large zt/U, the ratio t to V is the only important factor 
determining the phase (SF or SS, not the PS), because 
the lattice does not play an important role except at half 
or unit filling, and N affects only the SF-SS phase bound- 
ary (V affects the phase more strongly for larger N) . As 
a result, no rapid change in the ground state properties 
occurs around zV/U = 1 in contrast to the case of small 
or intermediate zt/U. 

The dot-dashed curve of perturbation 1 does not agree 
with the full numerical calculation (solid curve) because 
perturbation 1 exhibits the particle-hole symmetry of the 
hard-core model, which is distinctly broken here. How- 
ever, the dot-dashed curve of perturbation 2 is in ex- 
cellent agreement with the solid curve except at large 

Figure 12 shows the V dependence of Sn in the SS 
phase of the soft-core model for zt/U = 1. Sn is larger 
for larger TV at the same zV/U because the critical value 
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(a) 



1.10 



1.09 




v — perturbation 3 \ 

\ [SF-PS(SF+S 2 )] \ 
\ 



0.9 



0.95 
N 



(b) 



FIG. 8: Phase diagrams of the soft-core model for zt/U = 
0.3. (a) shows Entire phase diagram. Solid curve shows the 
continuous SF-PS(SF + Si) and SF-SS transitions. Dot- 
dashed curves are the results of the perturbative calculations 
Eq. [9] (perturbation 1), Eq. IB16I (perturbation 2), and Eq. 
IB211 (perturbation 3) for the SF-PS(SF + Si) transition. 
Two-dot-dashed lines at half filling and unit filling represent 
the solid phases (Si and S2, respectively), (b) Expansion 
of the region around zV/U = 1 and N = 1. Solid 
curve shows the continuous SF-PS(SF + S2) phase transi- 
tion; dashed curve shows the discontinuous SF-SS and PS(SF 
+ S2)-SS phase transitions. Dot-dashed curve shows Eq. [B24l 
(perturbation 3) for the SF-PS(SF + S 2 ) transition. 



of zV/U for the SF-SS transition is smaller for larger N 
(Fig. 11). Unlike the cases of zt/U = 0.03 (Fig. 6) 
or zt/U = 0.3 (Fig. 10), 6n is a smooth increasing 
function of N and zV/U, and the two curves of 6n for 
two different ./V values no longer intersect. Note that Sn 
does not change rapidly around zV/U ~ 1 as expected 
from the phase diagram (Fig. 11). 




1 1.2 1.4 1.6 

zV/U 



FIG. 9: V dependences of Sn = \Na—Nb | (solid curves) and 
ys 1 (dashed curves) in the soft-core model for zt/U = 0.3 
below half-filling. 




zV/U 



FIG. 10: V dependence of 8n — \Na — Nb \ of the soft-core 
model for zt/U = 0.3 above half filling. 



V. EFFECT OF IMPROVED CALCULATION 
ON THE ENERGY OF THE SOLID AND MI 
PHASES 

The energy of the solid and MI phases employed in 
the previous sections may be significantly higher than 
the exact energy. Hence, in this section, we improve the 
calculated energies in these phases by employing the per- 
turbation theory Q up to the order of t 2 /U or t 2 /V as 



Es 2 



zV 



2zt 2 
' U -V 
zt 2 



u 



2(z-l)V 
zt 2 



[2z-l)V-U' 



(12) 
(13) 

(14) 
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N 

FIG. 11: Phase diagram of the soft-core model for zt/ U = 1. 
Solid curve represents the SF-SS phase boundary at which the 
phase transition is continuous. Dot-dashed curves show the 
results of perturbation 1 (Eq. O and 2 (Eq. IB16|) . Two-dot- 
dashed lines at half filling and unit filling show the Si and S2 
phases, respectively. 
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zV/U 

FIG. 12: V dependence of Sn = \Na — Nb\ in the soft-core 
model at N = 0.3, 0.5, 0,7, and 0.9 for zt/U = 1. 



We can obtain the phase diagram by employing these 
improved energies. Hereafter, we assume z — 6 for 
the 3D cubic lattice. However, note that this improve- 
ment may be excessively favorable for the solid and MI 
phases and unfavorable for the SF and SS phases. There- 
fore, a correction on the order of t 2 /U or t 2 /V should 
also be applied to the energies of the SF and SS phases 
close to the SF (or SS)-solid (or MI) phase boundary. 
In addition, the denominators of the equations describ- 
ing Emi, Es lf and Es 2 diverge for V = U, V = 
and (2z — 1)V = U, respectively. Hence, we exclude 
the Si phase around V = and the S 2 phase around 
zV/U = z/(2z - 1) = 6/11. We also neglect the PS 
into two solids (Si and S2) which is unlikely but indeed 




N 

FIG. 13: Phase diagram of the soft-core model for zt/U = 
0.03. Energies of the solid and MI phases are improved (see 
text). 




N 

FIG. 14: Phase diagram of the soft-core model for zt/U = 
0.3. Energies of the solid and MI phases are improved (see 
text). 



appears in this approximation. Figures 13 and 14 show 
the phase diagrams for zt/U — 0.03 and zt/U = 0.3, 
respectively. (The phase diagram for zt/U = 1 is 
the same as that in Fig. 11.) In both figures, the re- 
gions of separated phases become large, as expected. For 
zt/U = 0.03, the PS into the SF and MI phases [PS(SF 
+ MI) phase] appears. Furthermore, for zt/U = 0.3, 
the SS phase disappears below half filling because of ex- 
istence of the PS(SF + Si) phase. These results sug- 
gest that the solid and MI energies have somewhat over- 
improved. On the other hand, the regions of PS into the 
SS and solid phases become large above half filling, which 
resembles the results of a 2D QMC study by Sengupta 
and coworkers |31| . 
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VI. CONCLUSION 

In this work, we studied the hard-core and soft-core- 
extended Hubbard models by the Gutzwiller variational 
wave function. We adopted a canonical ensemble and 
a linear programming method to include PSs more di- 
rectly in our calculations. In the hard-core model away 
from half filling, we showed that the PS(SF + Si) and 
SS phases are degenerate above a critical value of the 
nearest-neighbor interaction V , which is consistent with 
a previous MF study [23| . 

Unlike the hard-core model, the soft-core model at half 
filling has a possible SS phase between the solid and SF 
phases and all the phase transitions are continuous. 

Away from half filling, the phase diagram depends 
drastically on the transfer integral t. For small zt/U, 
the shape of the SF region is similar to that of the hard- 
core model. The PS(SF + Si) phase appears below half 
filling and the SS phase appears only above half filling, 
as in the 2D QMC studies. For intermediate zt/U, the 
SS phase appears close to the PS below half filling as in 
the 3D QMC study. The phase diagram becomes simpler 
for large zt/U, where only the continuous SF-SS phase 
transition appears, and the critical value of zV/U at the 
phase boundary is a smooth decreasing function of N. 

The nearest-neighbor interaction dependence of Sn = 
\Na — Nb\, which shows the density wave order of the 
SS phase, is also interesting. For a small zt/U of 0.03, 
Sn is a discontinuous function of zV/U ; furthermore, the 
SS phase has a small Sn for small zV/U(< 1) and large 
5n(~ 2N) (N A ~ 2N and N B ~ 0) for large zV/U(> 
1). In addition, Sn is larger for smaller (larger) N for 
zV/U < 1 (zV/U > 1). For an intermediate zt/U of 
0.3, the behavior of Sn is similar to that for a small n 
zt/U — 0.03. In detail, however, unlike the case of 
zt/U = 0.03, the two curves of Sn for two different N 
values continuously intersect around zV/U ~ 1.1 because 
Sn is a continuous increasing function of zV/U except at 
large N. For a large zt/U = 1, Sn is a smooth increasing 
function of N and zV/U, and the two curves of Sn for 
two different ./V values no longer intersect. 

Throughout this paper, we found that our perturbative 
calculations determined the phase boundary curves very 
well except for large 1). 

We also studied the effects of the improved perturba- 
tive calculation on the energy of the solid and MI phases. 
The improvement enlarges the region of PS into the SS 
and solid phases above half filling, which resembles the 
results of the 2D QMC study. However, the improve- 
ment is excessively favorable for the solid and MI phases, 
resulting in an unusual enlargement of the PSs. There- 
fore, the energies of the SF and SS phases should also be 
improved in future work. 

Because the Gutzwiller approximation is not precise, 
exact numerical calculations such as QMC simulations 
are needed to check our results. Although some details 
(such as the very complicated phase diagrams for N ~ 1 
and small zt) might be artifacts of our approximation, we 



believe that the important results in the phase diagrams 
and Sn values in the SS phase are worth studying in de- 
tail. For instance, the transfer integral dependence of the 
entire phase diagram seems to remain an open question 
not only in the 3D case but also in the 2D case, especially 
for the large zt/U regime, in which an SS phase below or 
at half filling might exist. 
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Appendix A: PERTURBATIVE CALCULATIONS 
FOR THE HARD-CORE MODEL 

1. SF-SS TRANSITION 

In the Gutzwiller approximation, the energy expecta- 
tion value of the Hamiltonian (Eq. [1]) is obtained as 

E = E kin + E^ t , 

-Ekin = — ztCAQCBOCAlCBli 

El t = ^4i4i. (Al) 

Here, we assumed that ca(b)u i n — 0, 1) is real without 
the loss of generality and used = ca(b)oCa(b)i 

and {a\( B ) a A(B)) — c a(b)i- ^° consider a possible SS 
phase infinitesimally close to the SF-SS phase boundary, 
we set ca(b)u = c n + Sa n (S/3 n ), where Sa n and Sj3 n are 
infinitesimal quantities. The normalization condition of 
the wave function for sublattices A and B is written as 

(c + 5a Q ) 2 + (ci + Sai) 2 = 1, 
(co+<%) 2 + (ci+<%) 2 = 1. (A2) 

The boson number condition is written as 

(ci +S ai ) 2 + (ci +<5/?i) 2 = 2N. (A3) 

These equations can be rewritten using the normalization 
(Sn c n = 1) an d boson number conditions (J2 n nc n — 
N) in the SF phase as 

4:C Sx + Sx 2 , + Sy = 0, 
AdSxi + 8x\ + Sy\ = 0, 
c Sy + ciSyi + Sx a Sy = 0. (A4) 

Here we introduced new variables Sx n (5y n ) = 5a n + 
{—)S(3 n and 8y n ^ 0. These equations can be further 
rewritten in the lowest order as 

AcqSxq + Syl = 0, 
4ciSxi + 5y 1 = 0, 
2co%) + 2ci(fyi = 0. (A5) 
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Hence, the value of 8x n is of the same order as Sy 2 , and 
Sx n can be rewritten in terms of <5y 2 . The energy expec- 
tation values can also be rewritten as 



kin 



-Zt 

zv r 

~2 



1 



1 



N(l-N)--(1-N)8yi--N5yi 



By using Eq. [A4 



TV 2 - Acfy 2 



(A6) 



zt 



E = -zt N(l - N) + yiV 2 



N 2 + (l-N) 2 
1-N 



zVN 



Sy 2 . (A7) 



If the coefficient of 5y\ is positive (negative), the phase 
is the SF (SS). By setting the coefficient of Sy\ equal 
to zero, we obtain the critical value of V for the SF-SS 
transition 



Vc 



N 2 + (l~N) 2 



N(l - N) 



(A8) 



This is nothing but Eq. |9] 



zV/2). To obtain the phase boundary, we substitute Eq. 
EU] into Eq. [ATT] The result is 



E 



-zt N(l-N) + —N 2 



-N + E< 



SM 



+ {ztN 2 + Z -^i 

+ [- zVN + zt(l - 2N)]N SM } 1SM (A12) 

in the lowest order of 7sm- For the PS that splits into 
the SF and Si phases, N$m = 1/2, and E$m = 0. 
If the coefficient of 7sm is positive (negative) , the phase 
is the SF (PS). By setting the coefficient of 7sm to zero, 
we obtain the critical value of V for the PS. The result 
is the same as Eq. |H1 In contrast, for the PS that splits 
into the SF and MI phases, N SM = 1, E SM = zV/2, 
and the coefficient of 7sm is positive definite: 

(W^)(A-1) 2 >0. (A13) 
Hence, the PS into the SF and MI phases does not occur. 



Appendix B: PERTURB ATIVE CALCULATIONS 
FOR THE SOFT-CORE MODEL 



2. PHASE SEPARATION INTO SF AND SOLID 
OR MOTT PHASE 

To study the PS from the SF phase into the SF phase 
and the Si (MI) phase, we set the ratio of the Si (MI) 
phase as 7sm- When the total boson density is N, the 
number density condition is written as 



(1-7sm)(A''-<W)+7sm-/Vs M - N. 



(A9) 



Here A^ SM = ^V Sl - 1/2 (A SM = A^mi = 1). If 
5N = 7sm = the phase is the uniform SF phase. In 
contrast, if 5N, jsm ^ 0, the phase is the PS consisting 
of the SF and Si (MI) phases, and N—8N is the number 
density of the SF phase. From Eq. IA9[ we obtain 



SN ~ (ATsm-AOt! 



SM 



(A10) 



near the SF-PS phase boundary (SN <C N and 7sm <C 
1). Because | ci| 2 = N — SN and |c | 2 = l-|ci| 2 , the 
system energy per site is 

E = 7sm-E : sm + (1 — 7sm)-Esf, 
E SF = -^|c | 2 |ci| 2 + ^|ci| 4 



-zt(l - N + SN)(N - SN) 



zV 



(N - SN) 2 



(All) 



Here, Eg F is the energy of the SF phase and Eg M is that 
of the Si (MI) phase: E SM = E Sl = (E SM = E M i = 



1. SS-Si TRANSITION AT HALF FILLING 

To obtain the SS-Si phase boundary at half filling (Eq. 
ITU]) , we set c^o = Sa , cai = 1 — Sax, ca2 = Sa 2 , 
cbo = 1 — 8f3o, cbi — SP\, and cb2 — 5^2 infinitesi- 
mally close to the phase boundary. If Sa n — 5f3 n = 
(n = 0, 1, 2), the phase is SF, but if 5a n ,5f3 n ^ 0, the 
phase is SS. The normalization conditions of the wave 
function for sublattices A and B lead to the following 
equations in the lowest order: 



Sax 



(Sal + Sa l)> 



SPl) 



(Bl) 



In addition, the boson number condition leads to 



Sax = Sa 2 + ^S(3'f + Sf3 2 . (B2) 



From these equations, we can eliminate Sao, Sax, and 
5/3q and obtain the energy in terms of Sa 2 , 5 fix, and 502 
as 



E 



U . 2 zV 

—Sai H 

2 2 2 

— zt ( \f2Sa2 



sp 2 



zV)S/3 : 



U 



zV 



-^ + -+{- + zV)y< 



U 



-zt 



(V2x+ 



+ S/3 2 +2Sf3 2 )S(3 



2y 



(B3) 
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where we have set #a 2 = x5/3i and <5/3 2 = y&fii- The 
minimization conditions of Eq. IB3I are written as 



ztx 



1 9E TJ r— 

7^2^- = Ux-V2zt 

8 Pi dx ^l + x 2 + 2y 2 



0,(B4) 



wir - <p+**y)v- h 2 ? 92 =°( B5 ) 



The latter equation leads to y = 0, because we assume 
small zt(<C ?/). From the former equation, we have 



y/2zt 

u 



zt 
U 



O 



zH 3 
U 3 



(B6) 



where 



X 



«o 



»02 



d2 = 



(V2C2 



■ c 



V c 



-2c 2 (V2c 2 + c ) 



( 2 + ^) C ^ + t^ C2 + C °) 



24 

.2 
1 



72c? + ^(V2c 2 



2 

2(\ : 



c 



-cf(l 



-co; 

V2C2 
CO 



cf 



(B9) 
(BIO) 
(Bll) 

(B12) 



By substituting the a; value obtained above into Eq. IB31 
we obtain the optimized coefficient of S(3f . If it is positive 
(negative), the phase is Si (SS). By setting it to zero, 
we obtain the critical value of zV for the SS-Si phase 
transition (Eq. [T0|) . 



2. SF-SS TRANSITION 

As explained in the text, we limit the Hilbert space to 
three states |0), |1), and |2). Furthermore, we set CAn = 
c n +6an and cb„ = c n +5(3 n where 5a n and S/3 n (n = 0, 
1, 2) are infinitesimal quantities. We determine the c n 
(n = 0, 1, 2) value required to minimize the energy of the 
SF phase (Sa n — 5f3 n — 0). As in the hard-core model 
(Appendix A.l), we introduce 6x n (6y n ) — 5a n +(—)5f3 n , 
where Sy n ^ for the SS phase. The wave function 
normalization and boson number conditions lead to 



cqSxq + aSxi + c 2 Sx 2 = ~^(Syl + Sy 2 + Sy 2 ), 



Sy 2 

2ciSx 1 + 4c 2 5x 2 = — — - Sy 2 , 
co^yo + c 1 Sy 1 + c 2 6y 2 = 



(B7) 



in the lowest order of Sx n and Sy n (Sx n are of the order 
of 5y 2 ). By eliminating Sxq, Sx±, and 8y1 in the above 
equations, we can write the energy per site in terms of 
Sx 2 , Syo, and Sy 2 as 



E = 

Ekin — 



E w 



Ekin 
-Zt 

1 



E; r 



c 



XSx 2 



(aady 2 + ao 2 <5z/ <5y 2 + a 2 8y 



E> 



II 

zV 
~2~ 



+ C 2 SX2 + ■ 

(cl+2cl) 2 



Sy 2 2 



(c Sy Q - c 2 Sy 2 



(B8) 



The total coefficient of Sx 2 for E in Eq. IB8I is zero, be- 
cause the variation in 8x n corresponds to that in c n in 
the SF phase and the values of c n were already deter- 
mined to minimize E. To be concrete, the energy of the 
SF phase is represented by c 2 when Sy n = as 



E SF = -zt{N -2c 2 )[V2c2 + ^l-N - 
+ 7 ^N 2 + Uc 2 2 , 



(B13) 



where Co and c\ were already eliminated by the the 
wave function normalization and boson number condi- 
tions. The minimization condition dEsp/dc 2 = leads 
to 



Uc 2 



zt 



c 2 



y/2, 



C2 



TV — 2c 2 
y/l - N + 6 



1-N + c 
v / 2c 2 



+ \ 1 -N + ct 



x \ c 2 + \/2(l- N + c 2 ) 



(B14) 



and we can easily verify that the total coefficient of 8x2 
for E in Eq. IB8I is zero. E is further rewritten as 

E = const. + A8yl + B5y 5y 2 + C5y%, 



A 
B 
C 



Zt zV 2 

J«o - -yCo, 

zt 

— a 02 + zVc c 2 , 

zt zV 2 U 
-a 2 



(B15) 



4 2 z 4 

Hence, Sy n becomes finite and the phase becomes SS 



when B 2 > 4AC, which corresponds to V > v£ F ~ 



where 



V, 



SF-SS 



— — 7—5 5 — : — — otrSlu J 

8 zt^ao + c c 2 a 02 + Cga 2 ) + UCq 



To obtain the numerical value of Eq. IB16[ we numerically 
calculate the value of c 2 that minimizes E (Eq. IB13[) 
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and determine Co by using the normalization and boson 
number conditions. By replacing c and c 2 in Eq. IB161 
we obtain the critical value V' ( ? F_SS . 

On the other hand, at half filling, we analytically ob- 
tain a power-law expansion of zt/U for Vq F_ (Eq. [TTj) . 
We start from the minimization condition obtained from 
Eq. IB13I up to the second order of c 2 : 

4^ = -zt(l + c, - 9c 2 ,) + 2Uc 2 = 0. (B17) 

dC2 

From the above equation, we obtain c 2 = zt/(2U) + 
z 2 t 2 1 \AU 2 ). We determine Co and c% by c 2 through the 
normalization and boson number conditions. Now, ao, 
ao2, and a 2 are written as 



From these equations and the minimization condition for 
Eqf (Eq. IB14|) . we determine the energy of the SF phase 
Esf in the separated phase as 



Esf = 



where 



-zt(N - 2c 2 2 )(y2c 2 

zV . 

T 



1 -N- 



+ ^N 2 + Ucj+YSN, 



Y 



-yc?(V2c 2 + c ; 

■ zVN ~T 



V2 

C2 



(B20) 



(B21) 



«2 



= 2 + 56 + US 2 , 
= \/2(l+2<f + 6<* 2 ), 

25 1 cx2 ^ 



= -d 

2 V 



(B18) 



respectively, where 6 = zt/(2U). By substituting ao, 
ao2, and a 2 into Eq. IB16[ we obtain Eq. 111! 



3. PHASE SEPARATION INTO THE SF AND 
SOLID OR MOTT PHASE 

As explained in the text and in Appendix B. 2, we limit 
the Hilbert space to three states: |0), |1), and |2). Wc 

assume that c„ = cl? (where n — 0, 1, 2) is determined 
to minimize the total energy of the SF phase (not the 
separated phase). The wave function normalization and 



boson number conditions are written as ^ ? 



.(0)2 _ 



and J2 r . 



,(0)2 _ 



= 1 

N, respectively. For the separated 

phase, we set Cj — > + Sc n and N — > N — SN for the SF 
phase in the separated phase and 7sm as the ratio of the 
solid (or MI) phase to the entire system. As in a similar 
calculation in Appendix A. 2, the boson number condition 
in the entire system is written as SN = 7sm(-/Vsm — 
N) in the lowest order of 7sm, where ./Vsm = 1/2 (1) 
for the solid Si (S 2 ) phase and ./Vsm = 1 for the MI 
phase. The normalization and boson number conditions 
for the SF phase in the separated phase are written as 



J2n(c { n ] +Sc n ) 2 = land£„n( C r+£c„) 2 = 
respectively. These equations are rewritten as 

Sc = -^-(2c 1 Sc 1 -SN), 
4c 

5c 2 = — l —{2 Cl 5 Cl +SN). 

4c 2 



(o) 



N-SN, 



(B19) 



Note that the 5c\ term is eliminated: the optimization 
process of Egp is similar to that in Appendix B.2, and 
the coefficient of 8c\ is found to be zero by Eq. IB14I 
Because the ratio of the SF phase in the separated phase 
is 1 - 7sm and SN — (N SM ~ N)~/ SM , the energy of the 
entire system E is written as 



E = £sf(1 - 7sm) + £sm7S 



M 



= -zt{N -2c 2 2 )[V2c 2 + sjl- N + 

+ ^-N 2 + Uc 2 
2 z 

+(Azt - BzV -CU + E SM ) 1SM , 



(B22) 



where 



.4 
B 
C 



2c 2 



= N 



(Nsm 



N 
~2 



„2 „2 

V2 C ° - ° 2 
coc 2 



(N SM - N) 



1 



= 4 + -(N SM -N), 



(B23) 



and Esm is the energy of the solid (MI) phase: E$m = 
{Esm — U/2) for the Si (S 2 ) phase, and Eqm = zV for 
the MI phase. If the coefficient of 7sm is positive (nega- 
tive), the phase is SF (PS). By setting the coefficient of 
7sm to zero, we find that the critical value of V for the 
SF-PS transition is 



Vi 



SF _ PS _ ztA ~ UC + Esm 
zB ' 



(B24) 
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